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ABSTRACT 

o 
o 

fS| . We present a Monte Carlo simulation for the response of the Laser Interferometer 

Space Antenna (LISA) to the galactic gravitational wave background. The simulated 
data streams are used to estimate the number and type of binary systems that will be 



< 



individually resolved in a 1-year power spectrum. We find that the background is highly 
non-Gaussian due to the presence of individual bright sources, but once these sources 



are identified and removed, the remaining signal is Gaussian. We also present a new 
estimate of the confusion noise caused by unresolved sources that improves on earlier 
C . estimates. 

o 

Subject headings: binaries: close — Galaxy: disk — gravitational waves — white dwarfs 
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^ 1. INTRODUCTION 

bJQ. 

Binary star systems are excellent sources of gravitational waves, and with roughly two-thirds 
of the ~10 11 stars in our galaxy in binary systems, there will be no shortage of targets for the 
proposed Laser Interferometer Space Antenna (LISA) (Bender et al. 1998). Binaries with periods 
less than a day may potentially dominate the response of the LISA observatory. Indeed, it is likely 
that the main source of noise for LISA over a portion of its band will be unresolved gravitational 
wave signals from galactic and extra-galactic binary star systems. Several studies (Evans et al. 
1987; Hils et al. 1990; Bender & Hils 1997; Postnov & Prokhorov 1998; Nelemans et al. 2001) have 
sought to model these populations and estimate their contribution to the gravitational wave power 
spectrum. In some instances these estimates have been combined with data analysis considerations 
to make predictions of the confusion noise caused by unresolved sources. Although a consensus has 
not formed on the expected background level, it is generally accepted that a galactic gravitational 
wave background does exist inside the LISA band. The difficultly in developing strong limits on the 
background level originate from the expectation that the background will be dictated by compact 
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binaries, that is, binary systems that contain white dwarfs, neutron stars, and stellar mass black 
holes. Since the electromagnetic luminosity of compact binaries is low, not enough sources have 
been observed yet to build reliable models for the populations. 

If current estimates of compact binary populations reasonably represent the true nature of 
the galaxy, then the superposition of gravitational wave signals from these populations will form 
a confusion limited background in the LISA band. That is, there will be enough sources that the 
received signals will interfere with each other to the point where individual binaries cannot be 
resolved (Crowder & Cornish 2004). When this occurs the background becomes a source of noise in 
the detector. However, around ten thousand systems will be resolvable due to either their isolation 
in frequency space (for sources with frequencies above ~3 mHz) or their relative brightness (for 
sources below ~3 mHz). 

Here we present a Monte Carlo model for the galactic gravitational wave background. Our 
goal is to better understand the role played by the rare, bright sources that dominate the observed 
signal, and to provide a more realistic level of the confusion background due to unresolved compact 
binaries. Our investigation of the galactic gravitational wave background is done in two phases. 
The first phase is to build a Monte Carlo simulation of the background by modeling each binary 
and processing the corresponding gravitational wave signal through a model of LISA. To reasonably 
represent the individual types of binaries we follow the population models presented in Hils et al. 
(1990) and Nelemans et al. (2004), which from here on will be referred to as HBW and NYZ 
respectively. Though less up to date, the HBW model has the advantage of being expressed in 
terms of explicit distribution functions, which allows us to generate multiple realizations. The 
more modern NYZ model employs a population synthesis code that has not been made public, so 
we were not able to generate our own realizations. Gils Nelemans was kind enough to send us 
a realization of the NYZ model to work with. The second phase of our study is to statistically 
characterize features of the galactic background as they are observed by LISA. 

As part of the modeling phase, the signal from each source is run through a realistic model 
of the LISA instrument response to produce simulated interferometry data. In doing so, it was 
necessary to develop a new, fast algorithm for computing the detector response in order to process 
the ~10 8 sources modeled in each realization. The algorithm is based on the frequency domain 
approach developed by Cornish & Larson (2003), and extended to include the detector transfer 
functions, arbitrary observation times, and frequency evolution of the sources. For the current 
study we set the observation time at T b s = 1 year. Examples of our simulated LISA data streams 
can be found at the Mock LISA Data Archive 1 . 

Using the simulated detector data we investigate issues that are of importance to the de- 
velopment of future data analysis algorithms. Among the quantities we investigate are tests of 
Gaussianity in the distribution of Fourier coefficients before and after bright sources are removed, 
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the number and type of bright sources, and the density (in frequency space) of bright sources. Our 
interest in bright systems stems from the idea that they will be identifiable in the data streams, 
and thus removable. They will also be instrumental in using the real gravitational wave data to 
study galactic populations and galaxy evolution. 

To model the removal of bright systems we use an iterative procedure using a running median of 
co-added instrument noise and galactic signals as the effective noise level. Sources were considered 
bright if they had a signal noise ratio (SNR) greater than some threshold with respect to the effective 
noise level. We considered both optimistic (SNR = 5) and conservative (SNR = 10) thresholds. 
As the bright sources are regressed from the data, the effective noise level drops, allowing more 
sources to be resolved. After several iterations we are left with a residual signal that is our estimate 
of the galactic confusion noise. Previous estimates of the confusion noise were derived by setting 
a maximum source density, with the reasoning that it would be impossible to resolve individual 
sources when the number of sources per 1/T b s frequency bin exceeded some threshold. Here we 
took a different approach that is based primarily on SNR thresholds, but we also studied the effect 
of a source density cut-off. Rather than working with the total source density we based our cut-off 
on the density of resolved sources. In other words, we considered the possibility that there will be 
a maximum number of sources that can be resolved per frequency bin. We studied the effect of 
a resolved source density cut-off by performing the iterative removals with and without a cut-off 
on the number of sources that could be resolved per frequency bin (we set a limit of one source 
per four bins). For some models the source density cut-off had a significant impact, but for other 
models the cut-off made very little difference. Our estimate of the confusion level for the HBW 
model differs from that of Bender & Hils (1997) despite the fact that we use the same galactic 
model. Our estimate of the confusion level for the NYZ model agrees fairly well with the Barack &; 
Cutler (2004) estimate. In both cases, our estimate is lower at low frequencies (below 1 or 2 mHz 
respectively), and higher at high frequencies. The differences are due to our differing approaches to 
modeling the signal identification and regression. We feel that our approach yields more realistic 
estimates. All our examples are for one year of observations, so the frequency bins have width 
A/ = 3.17 x 10~ 8 Hz. The level of the confusion background will drop for longer observation times 
as the sidebands get better resolved and the SNR increases. The reduction in the confussion noise 
over time means that fits to transient sources that occur in the first year of operation, such as 
a supermassive black hole merger, will continue to get better with time even though the source 
disappeared years ago! 

Recently Benacquista et al. (2004) and Edlund et al. (2005) simulated LISA time series for a 
population of galactic white dwarf binaries. While comparable to our approach, neither simulation 
was used to study source identification and subtraction. The statistical analysis of the background 
given in Edlund et al. (2005) focuses on the cyclostationary nature of the signal, whereas our 
statistical analysis focuses on simulating data analysis in order to better understand the galactic 
gravitational wave background. 

Since the study of the galactic gravitational wave background naturally divides itself into two 
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sections, modeling and characterization, the outline of the paper follows suite. Sections 2 and 3 are 
devoted to a description of the Monte Carlo simulation of the galactic close binary populations. It is 
here that we describe how the individual sources are modeled and convolved with a LISA response 
model. The next three sections calculate a number of statistical properties associated with the 
galactic background. Section 4 demonstrates that the galactic background is non-Gaussian in 
nature. In Section 5 we present our estimate of the confusion limited background and compare 
it to prior estimates. Section 6 describes the characteristics of the systems that are labeled as 
bright. The paper concludes in Section 7 with a discussion of the various assumptions used in the 
simulation and how changes in these assumptions may alter our results. 



2. GALACTIC MODEL 

The first step in building a Monte Carlo realization of the galaxy is to model an individual 
binary system. In general, a gravitational wave traveling in the k direction can be decomposed into 
two polarizations states, 

h(ct — k ■ x) = h + (ct — k ■ x)e + + h x (ct — k ■ x)e x , (1) 

where e + ' x are basis tensors used to describe the radiation's orientation. The scalar coefficients 
are referred to as the gravitational waveforms. For a circular binary with instantaneous angular 
orbital frequency f2, and component mass M\ and M<i the waveforms measured at the barycenter 
of the Solar System are 

h+(t) = A + cos(2V>) cos(20t + ip ) + A x sin(2i/>) sin(2fW + ip ) (2a) 
h x (t) = -A + sin(2V>) cos(20t + tp ) + A x cos(2^) sin(2fW + ip D ) , (2b) 

where the polarization amplitudes are given by 

2G 2 M 1 M 2 ( n 2 \ 1/3 , 2/ NN 

4G 2 M 1 M 2 ( ft 2 n 1/3 



= gr^ Uw+M,) ; cos(l> ' (3b) 

The angles ip and i describe the orientation of the binary as viewed by an observer in the barycenter 
frame, while ip is the initial phase. 

Gravitational waves carry away energy and angular momentum from the emitting system. 
Consequently a binary will slowly inspiral over time. For stellar mass galactic sources in the LISA 
band the period evolution can be adequately described by 

3/8 
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where M = (MiM 2 ) 3 / 5 (Mi + M 2 )~ 1 ^ is the so-called chirp mass and P Q is the initial orbital 
period. Equation (4) does make the assumption, which is used throughout this paper, that no 
other processes (e.g. mass transfer) are evolved in the binary evolution besides gravitational wave 
emission. 

Equations (2), (3), and (4) indicate that a circular binary is uniquely determined by a set of 
nine parameters: the component masses (Mi,M 2 ), initial orbital period (P ), binary orientation 
(ip,i), initial phase ((p ), and the distance to the source (r). Additionally, two angular variables 
(9, (f>) are used to locate the source on the celestial sphere. To model an individual binary requires 
an accurate representation of these nine parameters. 

The list of source parameters are separable into those that are extrinsic and intrinsic to the 
system. The extrinsic variables {r, 9, (ft, if), i, ip } do not influence the evolution of the binary. Instead 
they depend on the time of observation and on the location of the observer with respect to the 
binary. The remaining variables {M±, M2, P } directly effect the binary evolution through the 
emission of gravitational waves via equation (4). 



2.1. Extrinsic Parameters 

For the set of extrinsic variables there is a further separation into those that locate the source 
{r, 9, cf)} and those that describe the time of observation and orientation as viewed by a particular 
observer <p }. To derive a unique location for each source we use a cylindrically symmetric 
disk model of the galaxy with an exponential falloff in both the radial and vertical directions, 

p = Po e- r / r °e-\ z \ /z ° . (5) 

Here p Q is the space density at the galactic center, r a is the radial scale length, and z Q is the vertical 
scale height of the galactic disk. The values of r Q and z Q vary with the different types of binaries 
(i.e. cataclysmic variables, white dwarf binaries, etc.), but all types are assumed to obey the 
above model. The binary positions are simply described in galactocentric-cylindrical coordinates. 
The natural coordinate system for the LISA mission is heliocentric-ecliptic coordinates. Therefore, 
once the positions for the binaries are selected using the galactic position distributions, they are 
translated to the LISA coordinate system through a series of standard coordinate transformations. 

The observed orientation of a binary system is set by the principal polarization angle if) and 
the inclination angle 1. The inclination angle is defined as the angle between the line of sight to 
the binary h and the angular momentum vector of the binary L. The inner product of n and the 
angular momentum directions n ■ L is taken to be uniformly distributed between —1 and 1. The 
principle polarization angle describes the orientation of the semi-major axis of the projected binary 
orbit on the celestial sphere and is uniformly distributed between and ir. The distribution for 
(fo, which describes the positions of the binary components at time t = 0, is uniformly distributed 
between and 2ir. 
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2.2. Intrinsic Parameters 

The distributions for each intrinsic parameter {M±, M2, P } depends on the binary type under 
consideration. To model each of these parameters we used the distributions given in HBW. For 
this reason our galactic backgrounds include W UMa (3 x 10 7 ), cataclysmic variables (1.8 x 10 6 ), 
neutron star - neutron star (10 6 ), black hole - neutron star (5.5 x 10 5 ), and close white dwarf 
(3 x 10 6 ) binaries. The quantities in the parentheses indicates how many systems of that type are 
included in the simulation. Note that for most of our analyses we use the 10% reduced population 
of close white dwarf binaries as described in HBW as this allows us to to compare directly with 
prior results. 

We have elected not to include the unevolved binaries. The reason for this is that they are 
predominately very low frequency sources (< 1CT 5 Hz). As a result, their signals will be buried in 
the instrumental noise and, therefore, will not contribute to the observed galactic background. 



2.3. Barycentric Background 

Many of the prior studies of the galactic background approached the problem by estimating 
the net gravitational wave luminosity as a function of frequency in the Solar System barycenter. 
From the luminosity they then derive a gravitational wave strain amplitude using (Douglass k, 
Braginsky 1979) 

, (l6nG L gw \ 1/2 

where u> gw and L gw are the gravitational wave angular frequency and luminosity respectively. 

In order to make comparisons between our results and prior studies, we must relate the above 
expression for the strain amplitude to quantities calculated in our Monte Carlo simulation. To this 
end, we first note the relationship between gravitational wave luminosity and flux (Press & Thorne 
1972), 

L ""' - F »* = £n(>>l + K)> ( 7 ) 



4vrr 2 yw 16vrG 

where the angle brackets denote an average over several gravitational wave periods. Using this 
relationship the strain amplitude is rewritten as 

( 1 \ 1/2 
h=(-^(hl + hl)) . (8) 

From the waveforms given in equation (2) it follows that 

/I \l/2 

h=^-(Al + Al)j . (9) 
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The polarization amplitudes, A + and A x , are functions of the binary masses, distance to the source, 
orbital period, and inclination angle (see eq. [3]). 

Equation (9) gives the strain amplitude for a single source. To mimic a power spectrum in 
the Solar System barycenter frame we first bin the sources according to their frequencies. The bin 
widths are A/ = 1/T, where T is the total observation time (for our simulations T is set to one 
year). Once the sources are sorted the net strain amplitude per frequency bin is calculated from 



where Nb is the number of sources in the bin. Note that equation (10) accounts for constructive and 
destructive interference in the same way as in standard data analysis uncorrelated, random errors 
add quadratically with the square root taken after all errors have been included. Similarly for the 
background, the net strain amplitude is the square root of all individual polarization amplitudes 
added quadratically. 

Figure 1 compares our Monte Carlo results to HBW. (Note that the plots show spectral ampli- 
tudes hf, not strain amplitudes h. For monochromatic binaries the two are related by hf = \[Th 
where T is the observational period.) For each binary class we are in agreement. The smearing effect 
seen at high frequencies is due to empty bins in the spectrum. For the compact binaries the orbital 
period distributions have a small, but finite, probability at short periods. Consequently it takes a 
large number of draws against the period distribution to produce a source with an extremely short 
period. In the cases of the neutron star - neutron star and close white dwarf binaries the probability 
in the period distribution tails becomes small enough that for the number of sources included in 
our simulated background one would not expected a large number of extremely short period (high 
frequency) sources. This is why the HBW data extends beyond our simulated backgrounds. 

Figure 2 shows the total galactic background as viewed in the barycenter frame. The sharp 
rise in the background at / = 10 -46 Hz is due to the sudden increase in the number of W UMa 
binaries. The signal below / = 10~ 4 6 Hz is due to neutron star - neutron star binaries. If the 
unevolved binaries had been included, the galactic background would be approximately constant 
between 1 and 100 ^.Hz at a level of hf « 10~ 17 Hz -1 / 2 . 

As Nelemans et al. (2001) found with their population synthesis models of the galaxy, when 
the individual sources are modeled the background appears spiky. The large fluctuations in the 
background are due to a small number of bright sources. As we will show later, when these sources 
are removed from the background the spectrum becomes smooth. 

It is generally agreed that a confusion limited background arises when the average number of 
source per frequency bin is larger than unity. (Additionally the net strain per bin must also be 
larger than the intrinsic detector noise.) Figure 3 shows that the peak number of sources per bin 
is in excess of 10 5 near 0.04 mHz. The plot also demonstrates that for a large portion of LISA's 
spectrum the number of sources per bin is greater than ten. 
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Fig. 1. — A comparison of our full Monte Carlo simulation (gray) and a running average of the 
simulation (white) to those of HBW (squares). For each binary type we are in agreement. To ease 
comparison, the final figure shows a background using the full realization of white dwarf binaries. 
However, our simulated background uses the 10% reduce white dwarf population. 
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Fig. 2. — A realization of the galactic background as observed in the barycenter frame. The dark 
line is the all sky and polarization averaged LISA sensitivity curve (Larson et al. 2000). The 
jump at 10 -4 ' 6 Hz is due to sudden increase in the number of W UMa binaries. Had we included 
a realization of the unevolved binaries, the background levels would be roughly constant below 
10 -4 ' 6 Hz with a spectral amplitude of hj ~ 10~ 17 Hz -1 / 2 . 

3. DETECTOR BACKGROUND 

The backgrounds presented in the previous section were not convolved with a model for the 
instrument response. As with all fields of astronomy, the act of measuring incident radiation has to 
be properly understood in order to correctly interpret the signals. In this section we first summarize 
the LISA mission and then describe our models for the detector response. We then present the 
simulated galactic background as measured by the detector. 

LISA is a joint ESA/NASA mission planned for launch around 2015. The mission consists 
of three identical spacecraft in separate, slightly eccentric, heliocentric orbits inclined with respect 
to the ecliptic plane (Danzmann & Rudiger 2003). The orbits are carefully chosen such that the 
spacecraft constellation will form, and approximately maintain, an equilateral triangle with a mean 
spacecraft separation of 5 x 10 6 km. The center of the constellation, referred to as the guiding 
center, will have an orbital radius of 1 AU and trail the Earth by 20°. During the course of one 
orbit, the constellation will cartwheel once with a retrograde motion as seen by an observer at the 
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Fig. 3. — Number of sources per frequency bin. Again the dip in the number of sources at 10" 4 - 6 Hz 



is associated with the neglected unevolved binaries. 

Sun. LISA is sensitivity to gravitational radiation in the range of 1(T 5 to 1 Hz. 

The detector's motion introduces amplitude (AM), frequency (FM), and phase modulations 
(PM) into the gravitational wave signals (Cornish & Larson 2003). The amplitude modulation 
originates from the detector's motion sweeping the antenna pattern across the sky. The phase 
modulation results from the differing responses to each polarization state. The frequency (Doppler) 
modulation is due to the motion of the detector relative to the source. Since the bulk orbital and 
cartwheel motions both have a period of one year, the resulting modulations appear as sidebands 
in the power spectrum separated from the instantaneous carrier frequency by integer values of the 
modulation frequency, f m = 1/year. 

For our studies the detector response is modeled using a combination of the Extended Low 
Frequency Approximation, which is developed in the appendix, and the Rigid Adiabatic Approxi- 
mation as described in Rubbo et al. (2004). As explained in the appendix, to save on computational 
costs it is advantageous to simulate the detector response directly in the frequency domain where 
a quasi-monochromatic source will only have a small number of sidebands. At low frequencies, the 
bandwidth for a slowly evolving circular binary is 
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where R = 1 AU is the orbital radius of LISA and 9 is the colatitude of the source on the celestial 
sphere. For sources with gravitational wave frequencies below a few millihertz, the bandwidth is 
less than 100/ m and we can achieve a considerable saving in computational cost by working in 
the frequency domain. At higher frequencies the sources have larger bandwidths, more complex 
modulation patterns, and signals that evolve significantly in frequency, making them harder to 
model directly in the frequency domain. 

In the appendix it is shown that the Extended Low Frequency Approximation, which is ap- 
plied in the frequency domain, is only valid for frequencies below 7 mHz. For the few hundred 
signals with a carrier frequency above this cutoff, we use the more accurate time domain response 
model, the Rigid Adiabatic Approximation. To combine the results from each approximation, we 
first simulate the response for the signals that are above 7 mHz using the more detailed Rigid Adi- 
abatic Approximation and add them linearly in the time domain. We then perform a Fast Fourier 
Transform. For the sources that are simulated directly in the frequency domain, we coherently 
add the signals by summing the real and imaginary parts of their respective Fourier coefficients at 
each frequency. By adding the coefficients we maintain the phase information which dictates the 
constructive and destructive interference of the gravitational waves. The final detector response is 
the sum of the Fourier coefficients from the extended low frequency and adiabatic results. It is also 
at this time that we add in a detector noise realization using the prescription given in Rubbo et 
al. (2004). Figure 4 shows a particular realization for a Michelson signal in the frequency domain. 
Figure 5 shows the same signal in the time domain. 

In comparing the background as observed in the detector frame (fig. 4) versus in the barycentric 
frame (fig. 2) a striking feature is that it is lower by roughly a full decade across the entire spectrum. 
The reduction is due to two effects. The first effect is the detector efficiency, which relates the total 
signal power in the detector to the total signal power at the barycenter. The all-sky and polarization 
averaged detector efficiency is equal to -y/3/20 at low frequencies, and gets progressively worse at 
high frequencies. (This is immediately evident by comparing the sensitivity curve in figure 2 to the 
average Michelson noise in figure 4.) The second effect is due to the orbital motion of the detector. 
In the barycentric frame, most galactic binaries are well approximated as monochromatic. As LISA 
moves in its orbit the monochromatic signals arc modulated across multiple frequency bins. At 
high frequencies the spreading effects are evident by spectral power showing up in bins that were 
previously empty in the barycenter frame. At low frequencies the expectation is that the spreading 
from adjacent bins will cancel out due to the numerous sources found in each bin (see fig. 3). 
However, as will be shown in the next two sections, the galactic background is dominated by a few 
bright sources. When the bright sources are modulated there is not a compensating bright source 
in the adjacent bin. As a result the galactic background is also reduced at lower frequencies. 
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Fig. 4. — A realization of the HBW 10% galactic background as measured in the detector's frame. 
The dark line is the average Michelson noise associated with the detector. The galactic gravitational 
wave background is evident in the spiky structure between 0.1 and 10 mHz. 

4. STATISTICAL CHARACTER OF THE GALACTIC BACKGROUND 

Of great interest to the LISA mission is to broadly characterize the galactic gravitational 
background in a statistical sense. Such a characterization is essential to the development and 
implementation of data analysis algorithms which often make assumptions about the character of 
the noise. 

In the spectral regions of LISA's band where the galactic background dominates the detector 
response, the background becomes a source of noise. By inspection of the spectrum in figure 4, 
the galaxy is evident by the jaggedness between 0.1 and 10 mHz. Outside this region the galactic 
binary signals are weaker than the intrinsic detector noise. This is evident in the plot by the 
relative smoothness of the spectrum from bin-to-bin. One way to characterize the background is 
to statistically study the Fourier coefficient distributions in different regions of the spectrum. Of 
specific interest is finding out if the galactic background is characterized by a Gaussian distribution. 

Tests for Gaussianity are done using independent x 2 an d Kolmogorov-Smirnov tests. The 
Gauss tests are performed over a window of 512 bins and done at each frequency. Figure 6 shows 
the results of the Kolmogorov-Smirnov test performed on the real coefficients for the spectrum 
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Fig. 5. — A realization of the HBW 10% galactic background as measured in the detector's frame 
shown in the time domain (with no instrument noise). The annual modulation of the signal due to 
the detector motion is apparent by the two large lobes. 



shown in figure 4. The p values plotted along the ordinate axis are the probabilities that the set 
of measured deviates are Gaussian distributed. While it is impossible to say with certainty that a 
set of measured deviates are necessarily distributed with a specific distribution, they can be shown 
not to agree. Values for p above a hundredth are usually accepted as agreeing with the tested 
distribution function. For p values below a hundredth, which are seen in the millihertz region of 
the detector output, indicate that the measured set of data is not Gaussian distributed. 

Comparing the results of the Gaussian test (fig. 6) and the original LISA output (fig. 4) 
indicates that the detector response is non-Gaussian for frequencies in which the background is 
above the intrinsic detector noise level. Outside these regions, where the detector's noise dominates 
the output, the returned p values are consistent with a Gaussian distribution, as they should since 
the simulation of the noise is based on Gaussian distributions. Similar results are found for the 
imaginary coefficients, other channels of data, and when calculated using a \ 2 test- 
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Fig. 6. — The Kolmogorov-Smirnov Gaussianity test applied to the real Fourier coefficients for the 
detector output shown in figure 4. The presence of the galactic background is apparent by the low 
p values about / = 10 -3 Hz. Similar results are found for the imaginary coefficients, for the other 
data channels, and when using a \ 2 test. 

A common misconception is that the galactic background should be Gaussian distributed. 
This assumption is based on the Central Limit Theorem, which states that for a large sample of 
random deviates, regardless of their parent distribution, the distribution of average values will be 
approximately Gaussian. However, the Central Limit Theorem is not directly applicable to the 
galactic background since the net power in a single frequency bin may be dictated by a single 
bright source. Moreover, in a small region of the spectrum there are very few bright sources. Far 
less than the large number of data points required for the Central Limit Theorem to apply. 

5. CONFUSION LIMITED BACKGROUND 

5.1. Gaussian Nature of the Confusion Limited Background 

Large fluctuations originating from bright sources cause the tails of the expected Gaussian 
distributions to be enlarged. If the bright sources were removed from the data streams, then the 
remaining background would be Gaussian. To demonstrate this we first identify the bright sources, 
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subtract them, and then retest the remaining background for Gaussianity. 

Identifying all the bright sources in the actual LISA data streams is a difficult problem not 
yet fully solved. The modulation effects caused by LISA's orbital motion spread a source's spectral 
power across multiple frequency bins. Although the bandwidth over which a signal will spread 
is a known function of the gravitational wave frequency and sky position, if multiple signals are 
overlapping in a small region of frequency space, the true number of signals in the region may not 
be clearly identifiable. For our Monte Carlo models, when we generate each binary, the parameter 
values for each system are known. With this extra information we can quickly and accurately 
identify bright sources and remove them from the data streams. 

Our approach to identifying bright sources is to categorize them according to their signal-to- 
noise (SNR) ratio using the standard formula, 



where h(f) is the Fourier transform of the noiseless response to a single gravitational wave signal, 
and S n (f) is the one-sided noise power spectral density. A source is labeled as "bright" if its SNR 
is greater than 5 (optimistic) or 10 (conservative). The proper use of equation (12) requires a 
clear interpretation of the noise. We are interested in removing sources that are bright relative to 
the local power spectrum level. Therefore, the S n (f) curve must be a composite of the intrinsic 
detector noise and the galactic background. It is an effective noise for the detector, which we will 
emphasize by denoting it as S^f(f). 

To approximate the effective noise we calculate the median detector output. While representing 
the effective noise by the median response (as opposed to the mean response) is somewhat stable 
against bright sources, in regions near an extremely large SNR signal or where the density of bright 
sources is large, the median will still be influenced by these few signals. To account for the influence 
of the bright sources in calculating the (/) curve, we perform our calculations iteratively. We 
start with the median response of the initial output and calculate the number of bright sources 
with reference to this curve. We then remove the bright sources exactly using the same detector 
response approximation that generated them. The justification for using the exact parameters is 
a matter of simplicity, since data analysis algorithms are still being developed. Moreover, with a 
SNR threshold of 10 the errors in the recovered parameters will be small. The result of removing 
the bright sources is a new background from which we can calculate a new median response. From 
the new median response we calculate the number of bright sources with respect to the new (/) 
curve. We iterate this procedure until there are no new bright sources being identified. In some 
examples, such as with the 100% HBW model, the procedure does not appear to converge, so the 
subtraction was stopped after 10 iterations. 

Previous estimates (Bender & Hils 1997; Barack & Cutler 2004) of the confusion background 
ignored the relative brightness of the sources and focused instead on the source density. These 
estimates defined the confusion regime in terms of the number of sources per frequency bin. Outside 




(12) 
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of the confusion regime sources could be resolved and removed, while inside the confusion regime 
the sources acted as a source of noise. This notion of source confusion is based on linear algebra: 
The signal from each galactic binary is described by 7 or 8 parameters (Cornish &; Porter 2005), and 
there are 4 data points per frequency bin (two independent channels, each with a real and imaginary 
part). Thus, one needs at least 2 frequency bins per source to have as many data points as there are 
unknowns. While this estimate is very crude, it is unlikely that a data analysis algorithm can be 
found that beats the 0.5 source per bin saturation point by very much (methods such as Maximum 
Entropy (Jaynes 1995) introduce priors that can help tame under constrained systems, but they 
can only do so much). In order to study the effect of a source density cut-off we repeat the SNR 
based source subtraction procedure with a maximum resolved source density of 0.25 sources per 
bin using a 100 bin average. The density of one source per four bins was chosen as intermediate 
between the capabilities of existing algorithms (Cornish & Crowder 2005) and the saturation point 
described above. 

Table 1 lists the number of bright sources removed at each iteration and the total number 
of bright sources removed for several different realizations of the galactic background if no source 
density cut-off is applied. In each of the three HBW 10% realizations, the total number of sources 
was fixed at 3.635 x 10 7 . All three realizations gave similar results. The results of a HBW 
background using the full 100% of the white dwarves, along with the NYZ white dwarf background 
(2.6 x 10 7 sources) are also shown. For comparison we also include the optimistic case of using a 
subtraction threshold of SNR = 5. 

Table 2 lists the number of bright sources removed at each iteration and the total number of 
bright sources removed for several different realizations of the galactic background when a maximum 
density of 0.25 resolved sources per bin is applied. The HBW 10% and NYZ models are only slightly 
affected by the cut-off, while the HBW 100% and SNR > 5 version of the HBW 10% model are 
significantly affected. The impact of the cut-off is evident in Figure 7, where the effective noise 
levels for the HBW 100% (SNR > 10) and HBW 10% (SNR > 5) models are shown with and 
without the source density cut-off applied. 

Shown in figure 8 is the same Michelson signal as before (fig. 4), but with the the bright 
sources removed. Visual inspection of the spectrum shows that without the bright sources the bin- 
to-bin fluctuations are much smaller; an indication that the Fourier coefficients may be Gaussian 
distributed. Figure 9 confirms this hypothesis. When the bright sources are removed from the data 
streams the galactic background is Gaussian in nature. 

5.2. Confusion Limited Background Estimate 

The sources that are not flagged as bright will give rise to a confusion limited background that 
acts as an effective noise source for LISA. Taking the list of unresolved sources that remain after the 
simulated data analysis procedure described in the previous section, we can generate estimates of 
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Table 1. Number of New Bright Sources Identified at Each Iteration 



Iteration 


10%(1) 


10%(2) 


10%(3) 


100% 


10% (SNR=5) 


NYZ 


1 


6795 


6848 


6793 


10,736 


14,346 


8583 


2 


3806 


3723 


3712 


8084 


10,620 


4803 


3 


1693 


1803 


1846 


5340 


5480 


2090 


4 


817 


929 


866 


3550 


2577 


866 


5 


395 


457 


456 


2597 


1129 


473 


6 


226 


250 


225 


2014 


525 


266 


7 


114 


172 


123 


1628 


278 


191 


8 


59 


97 


76 


1400 


144 


125 


9 








1264 


72 




10 








1136 


38 




Total 


13,905 


14,279 


14,097 


37,749 


35,209 


17,396 



Table 2. Number of New Bright Sources Identified at Each Iteration with a Source Density 

Cut-Off Applied 



Iteration 


10% 


100% 


10% (SNR=5) 


NYZ 


1 


6795 


10,736 


14,346 


8583 


2 


3751 


7850 


4898 


4803 


3 


1669 


4440 


194 


2007 


4 


724 


1531 


20 


732 


5 


271 


463 


1 


325 


6 


79 


157 




186 


7 


30 


66 




79 


8 


12 


1 




48 


9 


7 






27 


Total 


13,338 


25,244 


19,459 


16,788 
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Fig. 7.— Effective noise levels for the the HBW 100% (SNR > 10) and HBW 10% (SNR > 5) 
models with (solid) and without (dotted) a source density cut-off of one resolved source per four 
frequency bins. 

the confusion noise in either the barycentric frame or in the instrument data channels. The former 
is useful for making comparisons with earlier work, while the latter is better suited to studying the 
effect of the confusion background on LISA's ability to resolve other types of gravitational wave 
signals. 

Figure 10 compares our barycenter and detector frame confusion noise estimates to the barycen- 
ter estimate of Bender Sz Hils (1997). We have multiplied our detector frame result by a factor of 
1^/20/3 to account for the average detector efficiency Our estimate is lower than the Bender & Hils 
(1997) estimate at low frequencies and higher at high frequencies. It is important to note that both 
estimates use exactly the same HBW model for the compact galactic binaries. The differences in 
the confusion noise estimates are due to the different way in which we modeled the data analysis 
procedure. Note that our results can only be compared below ~1.3 mHz. Above this frequency the 
Bender & Hils (1997) estimate is dominated by extragalactic sources, which we did not include in 
our simulation. 
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Fig. 8. — The same Michelson signal from before, but after the ~10 4 bright sources have been 
removed from the data streams. 



A simple piecewise fit to our confusion noise estimate of the HBW 10% background in the 
detector frame is given by 



'conf 



(/) 



' io- 


-45.9 J-2.6 


10" 


- 4A <f< io- 3 - 2 


10" 


-50.38 f-4.0 


io- 


- 3 - 2 < / < io- 2 - 8 


< 10" 


-78.38 J-14-0 


10" 


" 2 - 8 < / < 10- 2 - 65 m 2 Hz- 1 


10" 


-126.08 J-32.0 


10" 


-2.65 < j < IO -2 ' 55 


,10- 


-62.33 j-7.0 


10" 


-2.55 < f < 1Q -2.1 



Similarly, for the NYZ white dwarf binary background we found 



'conf 



(/) 



10" 


-44.62 j- 


-2.3 


10' 


" 4 < / < 10" 


-3.0 


io- 


-50.92 j-- 


-4.4 


10" 


3 < / < 10" 


-2.7 


10" 


-62.8 j- 


3.8 


10" 


- 2 - 7 < f < io- 


-2.4 


io- 


-89.68 f- 


-20.0 


10" 


- 2 - 4 < / < 10" 


-2.0 



m 2 Hz- 



(13) 
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These fits do not include instrument noise, and have been quoted in terms of position noise in order 
to avoid ambiguities in the path length scaling. (Our simulations normalize by the round trip path 
length, while other studies normalize by the detector arm length.) For comparison, our simulations 
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Fig. 9. — A running Gauss test applied to the detector output after the bright sources have been 
removed. Unlike before, all returned p values are consistent with a Gaussian distribution. 



of the Michelson response used an instrument noise spectral density 

1 



Sn(f) ~ ^ 2 



45 pos + 8(l + cos 2 (///,)) 



(15) 



(2vr/) 4 . 

with position noise S'pos = 4 x 10 -22 m 2 Hz" 1 and acceleration noise S acc \ = 9 x 10~ 30 m 2 s _4 Hz _1 . 
The choice of instrument noise level only has a weak effect on our results as the unresolved galactic 
background is the main source of noise from 0.1 mHz to roughly 3 mHz. 

A true confusion limited background is what remains after a full data analysis procedure has 
removed all identifiable signals. At present such an algorithm has not been fully implemented, 
though good candidates now exist (Cornish & Crowder 2005; Cornish et al. 2006). Our method 
of removal, by which we remove a source using the same response approximation that included it, 
mimics a true data analysis procedure, but it fails to include some of the subtle nuances associ- 
ated with signal and noise confusion (Cornish et al. 2006). Though we did not use a true source 
removal algorithm, we expect that figure 8 represents a good approximation to how the confusion 
limited background will appear in the LISA output. The best estimates of the confusion noise 
may ultimately come from Bayesian methods which treat the effective noise level as another model 
parameter to be estimated (Umstatter et al. 2005). 

Figure 11 shows the effective noise levels (confusion + instrument noise) for the HBW 10% 
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Fig. 10. — Our estimates for the HBW 10% barycenter (solid line) and detector frame (dot-dash 
line) confusion noise levels compared to the Bender & Hils (1997) estimate (dotted line). Our 
detector frame confusion noise has been multiplied by a factor of a/20/3 to account for the average 
detector efficiency. 

and the NYZ models. These estimates were produced using the conservative SNR = 10 criteria for 
bright source subtraction, and a maximum density of 0.25 resolved sources per bin. In contrast to 
what we found with the HBW 100% and HBW 10% (SNR > 5) models, the effective noise levels 
are little changed by the source density cut-off. Also shown is the effective noise estimate used by 
Barack and Cutler (Barack & Cutler 2004). The Barack-Cutler curve agrees quite well with our 
NYZ curve, though it does slightly overestimate the noise level between 0.5 and 2 mHz. 

6. BRIGHT SOURCE STATISTICS 



The bright sources represent signals that are identifiable in LISA's output. By understanding 
their location and separation (in the frequency domain) proper data analysis tools can be developed 
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Fig. 11. — Estimates of the effective noise level for the HBW 10% model (dotted line), and the 
NYZ model (solid line). The effective noise curve used by Barack-Cutler Curve (dashed line) is 
also shown. 

and applied in the search for their signals in the detector output. Also of interest are the properties 
of the bright sources. Are they near to us or distributed throughout the galaxy? Can we identify 
the type of binary system? The next two sections address these issues. 

6.1. Bright Source Density 

For issues concerning data analysis, an interesting quantity to know is the number of bright 
sources per frequency bin. Figure 12 is a plot of the average number of bright sources per frequency 
bin using a 100 bin window to calculate the average. The peak densities occur at ~2 mHz for the 
HBW 10% and NYZ models, whereas the HBW 100% model has a maximum at ~3 mHz. For 
the HBW 100% (SNR > 10) model and the HBW 10% (SNR > 5) model the maximum bright 
source density reached one source per bin. This is why the resolved source density cut-off of 0.25 
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Fig. 12. — The accumulated number of bright sources per frequency bin using a 100 bin average, 
both with (black) and without (grey) a bright source density cut-off. Starting in the upper left and 
going clockwise we have: HBW 10% with SNR > 10, HBW 10% with SNR > 5, HBW 100% with 
SNR > 10, and NYZ with SNR > 10. 

made such a big difference in those cases. The bandwidth of a typical source at this frequency is 
approximately twenty frequency bins (for one year of observation). As a result, in the peak density 
region there are bright sources whose power at least partially overlaps. 



6.2. Bright Source Characteristics 

An interesting question to ask is what property makes a particular binary bright? At high 
frequencies, where the number of sources per bin is less than unity, a source is considered bright 
if its signal is greater than the intrinsic detector noise. However, at low frequencies, where the 
number of sources per bin can be in the thousands (see fig. 3) , each source must compete against 
the other sources in the bin to become detectable. 
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log(//Hz) 

Fig. 13. — Distribution of outlier distances as a function of their frequency. Filled squares are 
white dwarf binaries, open circles are black hole - neutron star binaries, asterisks are are neutron 
star binaries, and open squares are cataclysmic binaries. 



To see what makes a source bright at low frequencies recall the functional form of the intrinsic 
amplitude, 

c 4 r \G(M 1+ M 2 )J ' 1 ° j 

Given a particular frequency bin the orbital frequency does not vary by more than a bin width, 
leaving only the masses and the distance to the source to determine if a signal is bright. 

Figure 13 shows the distance versus frequency for each bright binary coded by the type of 
system it is. Evident from this figure is the mass segregation associated with the bright systems. 
At low frequencies, where the number of sources per bin peaks, the more massive black hole - 
neutron stars are detectable. At higher frequencies the black hole - neutron star binaries are less 
numerous (see fig. 1) and the white dwarf binaries dominate the list of bright sources. 

To associate a binary to a particular type requires estimating the component masses. Un- 
fortunately for sources whose frequency evolution is too small to detect, there is a mass-distance 
degeneracy that prevents direct mass measurements (Cutler 1998). Table 3 shows the number of 
the bright sources that evolve by a measurable amount. The frequency evolution is considered 
measurable if the change in gravitational wave frequency during the time of observation is 
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greater than the width of a frequency bin (l/T Q b s ). 

It is interesting to note that very few sources will have a measurable frequency evolution during 
one year of observation. The two main reasons for this is that the frequency evolution is dictated by 
mass and initial frequency (see eq. [4]). The massive black hole - neutron star binaries are located 
at low frequencies, while the higher frequency white dwarf binaries have smaller chirp masses. 

Figure 13 also demonstrates that we should be able to see individual binaries distributed 
throughout the galaxy. This introduces the tantalizing prospect of using gravitational wave data to 
map galactic populations. However, since there is a mass-distance degeneracy for systems that do 
not evolve appreciable during the lifetime of the detector, only a small fraction of the identifiable 
sources can be used for such an analysis. 

7. SUMMARY AND CONCLUSIONS 

In this paper we have presented a Monte Carlo simulation of the galactic gravitational wave 
background as it would be detected by the proposed LISA mission after one year of operation. For 
the intrinsic binary properties we used the distributions given in Hils et al. (1990) and Nelemans et 
al. (2004). Our key findings are: the galactic background levels will be reduced in the detector frame 
as compared to the barycenter frame, prior to the removal of the bright sources the background is 
not characterized by a Gaussian distribution, and of the ~10 4 identifiable sources only ~10 2 are 
evolving and thus identifiable by type. We have also derived a new estimate for the confusion limited 
background that differs from other estimates given in the literature. Below 1 mHz our estimate is 
lower than previous estimates, while above 1 mHz our estimate is higher than previous findings. Our 
calculation of the bright source density (in frequency space) suggests that data analysis algorithms 
can be developed that are capable of resolving ~10,000 galactic binaries from a one year LISA data 
stream. Of these, roughly three hundred will be measurably chirping, allowing the determination 

Table 3. Number of Subtracted Binaries by Type 



Realization 1 Realization 2 Realization 3 



Type 


Included 


Removed 


Evolving 


Removed 


Evolving 


Removed 


Evolving 


W UMa 


3 x 10 7 




















CB 


1.8 x 10 6 


1 





2 





1 





NS - NS 


10 6 


11 





11 


2 


7 





BH - NS 


5.5 x 10 5 


591 


10 


585 


12 


550 


8 


WD - WD 


3 x 10 6 


13,302 


346 


13,681 


353 


13,539 


337 


Total 


3.635 x 10 7 


13,905 


356 


14,279 


367 


14,097 


345 
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of the chirp mass and the distance to these sources. The number of resolvable galactic sources, 
especially the number of measurably evolving systems, will increase significantly after several years 
of observation (the resolution of the frequency derivative improves with observation time, T G b s , as 

^obs )■ 

An important point to keep in mind about the results presented here is that they assumed 
particular descriptions for the galactic distribution, total number, and source characteristics of 
each binary population. A different collection of models of the extrinsic parameters may return a 
slightly different set of results as is seen in the NYZ backgrounds. However, the main conclusions 
drawn here are largely dictated by three quantities, the total number of binary systems, their period 
distributions, and the component masses. 

The number of binary sources in the galaxy can impart a noticeable difference in the back- 
ground levels. As radiation from the binaries converges on the detector, the random phase differ- 
ences will cause constructive and destructive interference. Statistically the problem is analogous 
to a random walk. As a result, the net spectral amplitude per frequency bin will grow as y/~N. 
For every factor of one hundred difference in the number of sources, the background levels raise 
or lower by a decade respectively. For electromagnetically visible binaries (W UMa, unevolved, 
and cataclysmic binaries) galactic surveys have placed stringent constraints on the total number of 
such binaries. However, for the compact binaries (neutron star - neutron star, black hole - neutron 
star, and white dwarf binaries) equivalent surveys have failed to place strict bounds on the total 
numbers. One return of the LISA mission will be to place limits on the populations by measuring 
the galactic background median levels and the number of bright sources from each population. 

While the overall level for the galactic background level is partially influenced by the total 
number of systems, the background level is also dictated by the component masses via equation (16). 
By comparing our HBW based simulations to the population synthesis approach of Nelemans 
et al. (2004), we have found that relatively small changes in the component masses can lead to 
significant changes in the background levels. Current uncertainties in the true mass distributions 
arise form a lack a observational data and a theoretical understanding of mass transferring stages 
during formation. The NYZ white dwarf binaries are typically composed of two light components, 
M < O.5M , whereas the HBW white dwarf binaries are typically composed of one light component 
and one heavy component, M > M Q . This results in the HBW systems having mean chirp masses, 
M = (M 1 M 2 ) 3/5 /(M 1 + M 2 ) 1/5 , almost a factor of two higher than the NYZ systems. Since the 
amplitude of signals scales as the difference in component masses translates into a factor of 

~ 3.2 increase in the HBW background relative to the NYZ background for the same source density. 
It is an interesting numerical coincidence that the 10% reduction in source density proposed by 
HBW yields a \/l0 ~ 3.2 reduction in the amplitude of the background, resulting in a background 
level comparable to the NYZ model. Despite this fortuitous cancellation, the HBW 10% and the 
NYZ models yield different confusion noise estimates due to the reduced number of sources in the 
HBW 10% model. 
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The period distributions will also impart a noticeable change in the background. If, for ex- 
ample, there is a mechanism that suppresses low period white dwarfs, the values in table 6.2 and 
figure 13 would change. Conversely, we can invert the problem and ask questions such as what 
would the high frequency end of the galactic background look like if certain physics are included in 
the models for close white dwarf binary production? 
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Agreement NCC5-579. L JR and SET acknowledge the support of the Center for Gravitational Wave 
Physics. The Center for Gravitational Wave Physics is supported by the NSF under Cooperative 
Agreement PHY 01-14375. SET is also supported by a Pennsylvania State University Graduate 
Fellowship. 



A. DETECTOR RESPONSE 

For a gravitational wave traveling in the k direction we can express a single channel LISA 
response as a linear sum of responses for each polarization state, 

s(t) = A + F + (t) cos $(t) + A X F X (t) sin $(t) , (Al) 

where the wave phase is given by 

$(t) = 2irf + irf t 2 + y> -<f> D (t). (A2) 

Here /„ is the initial gravitational wave frequency, f Q is the initial frequency derivative, and A+ x 
are the polarization amplitudes given in equation (3). The Doppler modulation of the signal is 
given by 

* D (t)2L^k-Xi(t). (A3) 

where we have neglected the small correction due to the frequency evolution f . The antenna beam 
pattern functions F +,x (t) describe LISA's time varying sensitivity to each polarization and are 
given by 

F + (t) = ^ ( cos(2</>)L>+(i) - S m(2ip)D x (t)) (A4a) 

F x (t) = i(sin(2V>)-£> + W + cos(2V>)£> x (t)) , (A4b) 

where the two-arm detector response functions are defined as 

D+(t) = d%%j{t, f gw ) - d+T tk (t, f gw ) (A5a) 
D x (t) = d x 3 % 3 (t^ gw ) - d x k T lk (t,f gw ) (A5b) 

(A5c) 
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with 

d±(t) = (f ij (t)®f ij (t)):e + (A6a) 
dijit) = (f ij (t)C^f ij (t)):e x . (A6b) 

The colon denotes a double contraction, x : y = x ab y a b, with repeated indices implying a sum- 
mation, fij(t) is a unit vector that points from spacecraft i to spacecraft j, and Tij(t,f gw ) is the 
round-trip transfer function for the arm connecting the i and j spacecraft, 



Tij(t,f) = - 



sinc {Jj- (l - fc • M*))) ex P ( _ *^T (3 + ■ M*))) 
+ sine f X (l + ■ hj{tj)\ exp (l + k ■ r^(t)) 



(A7) 

The quantity /* = c/2itL is referred to as the transfer frequency which for LISA (L = 5 x 10 6 km) 
has a value of 9.54 mHz. The transfer frequency is approximately the point at which a gravitational 
wave will "fit inside" the detector arms. 

The polarization basis tensors are expressed in terms of two orthonormal unit vectors, 

e + = u® u — v ® v (A8a) 
e x =u®v + v<g>u. (A8b) 

The unit vectors u and v, along with the propagation direction of the gravitational wave, k, form 
an orthonormal triad, which may be expressed as functions of the source location on the celestial 
sphere (9, <f>), 

u = cos(0) cos(0)x + cos(0) sin((j))y — sm(6)z (A9a) 
v = sin(0)x — cos(0)y (A9b) 
k = — sin(0) cos(0)x — sin(0) sm(4>)y — cos(9)z . (A9c) 



To calculate the unit vectors (t) we use the spacecraft coordinates given by Cornish & Rubbo 
(2003), 

x(t) = i?cos(a) + ^ei?(cos(2a - (3) - 3cos(/?)) (AlOa) 

y(t) = Rsm(a) + ^ei?(sin(2a - (3) - 3sin(/3)) (AlOb) 
z(t) = -V3eRcos(a - 0) . (AlOc) 

In the above R = 1 AU is the orbital radius of the guiding center, e is the orbital eccentricity, 
a = 2irf m t + k is the orbital phase of the guiding center, f m = 1/year is the modulation frequency, 
and (3 = 27r(n — l)/3 + A (n = 1, 2, 3) is the relative phase of the spacecraft within the constellation. 
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The parameters k and A give the initial ecliptic longitude and orientation of the constellation 
respectively. Note that to linear order in the eccentricity, which is the order we work to, the 
triangular formation is rigid with arm lengths given by L = 2\f3eR. 

Equation (Al), and the relationships that follow it, represent the Rigid Adiabatic Approxima- 
tion described in Rubbo et al. (2004). An important property of the approximation is that it is 
implemented in the time domain. To properly model sources with frequencies up to 10 mHz re- 
quires a minimum of ~10 6 data points for an observational period of one year. The computational 
cost of simulating the response for over 10 7 sources in the time domain is prohibitive. 

A desirable alternative is to work directly in the frequency domain. The advantage of doing 
so is that the number of relevant Fourier coefficients is small. Here the concept of "relevant" are 
those coefficients that contain a high percentage (~98%) of the spectral power. For reference, a 
moderate signal-to-noise ratio source at 10 mHz, observed for one year, will spread across sixty-five 
frequency bins. It is possible to derive an analytic expression for the Fourier transform of the time 
domain signal, equation (Al). The calculation has four steps. First, the amplitude and frequency 
modulations are decomposed into harmonics of the detector's orbital frequency f m . Second, the 
frequency evolution term, exp(7ri/ D t 2 ), is Fourier transformed. Third, the product of the orbital 
harmonics and the barycenter wave function exp(27rz/ t) are Fourier transformed, yielding a Fourier 
series who's coefficients are products of the harmonic amplitudes and the cardinal sine function. 
Finally, the Fourier series from steps two and three are convolved to give the complete finite time 
Fourier transform of the time domain signal. Our calculation generalizes the expression derived by 
Cornish & Larson (2003) to allow for arbitrary observation times, chirping sources, and the effect 
of the instrument transfer functions. 

The first step in the calculation relies on the fact that the functions F + (t), F x (t), and $£>(£) 
owe their time variation to the orbital motion of the detector. Thus we may decompose each of 
these functions into harmonics of the orbital frequency f m , 

F+(t) = Y,Pne 2ninfmt (Alia) 

n 

F x (t) = Y,c n e 2ninfmt (Allb) 

n 

e «*z>(t) = ^2d n e 2ninfmt . (Allc) 

n 

The coefficients d n are given by the Jacobi- Anger expansion, 

d n = J n (27rf(R/c) sin 0) e **r/2-*) , (A12) 

where J n is the Bessel function of the first kind of order n. Deriving expressions for p n and c n is 
complicated by the transfer functions that appear in equation (A5). While it is possible to perform 
the harmonic decomposition exactly using the Jacobi-Anger expansion, a simpler approach is to 
Taylor expand the transfer functions in powers of ///* then decompose each term into orbital 



-30- 



harmonics. For our current needs a second order expansion is sufficient, 
%jitJ) = l-i<2 + k-f ij (t)) (J^j 

~\ (4 + 3 (2 + k- h 3 it)) + (2 + k ■ r lj{ t)) 2 ) [£) 2 + O ( £) 3 . (A13) 

We can now re-express %j, dfj, and d*j in terms of orbital harmonics: 

4 



n=— 4 
4 



4 (*»/)= E * ¥ v,ne Mn,mt 

n=— 4 
4 

= E dX U,n e27rm/m '- (A14) 
n=— 4 

Convolving these expansions yields 

4 4 

P»=E E (dtj,lTij,m ~ ^%k,m) (A15a) 

4 m=— 4 
4 4 

c ™ = E E " ^,i%k, m ) , (A15b) 

l=-4m=-4 

where n = I + m. The range of the sums in the above harmonic decomposition can be traced 
to the functional form of the fij(t) unit vectors. From equation (A10) we see that the harmonic 
decomposition of rij{t) involves a sum from -2 to 2. The one-arm detector response functions and 
the expanded transfer functions are each a function of f 2 and, therefore, their decompositions range 
from -4 to 4. 

The second step is to Fourier transform the frequency evolution term, which yields Fourier 
coefficients q n that can be expressed in terms of error functions of a complex argument. The third 
step is to Fourier transform the remaining terms over the finite observation time T which has the 
effect of introducing the cardinal sine function, sinc(x) = sin(x)/x. Putting everything together we 
have 

s 3 = 2 e ^° E * E sinc(x im )e^™ £ [a +Pii + e^l 2 A x c n ) ^ d p , (A16) 

k l n p 

where xi m = Tr(lf m + foT — m), j = k + I and m = n + p. 

To test the range of validity for the new approximation, referred to as the Extended Low 
Frequency Approximation (ELF), we calculated a normalized correlation between the new approx- 
imation and the Rigid Adiabatic Approximation (RA), 

r(/ ) = {sRA(f)\s E LF(f)) (A17) 



V(4a(/))(4lf(/)) 
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Fig. 14. — The correlation between the Extended Low Frequency Approximation and the Rigid 
Adiabatic Approximation. 



The rigid adiabatic adequately describes LISA's response up to ~0.3 Hz (Rubbo et al. 2004) making 
it a valid benchmark to compare the new approximation against. Figure 14 shows the correlation 
between the two approximations for a randomly selected binary who's masses, sky location, etc., 
were held fixed while the correlation was tested at different frequencies. To safely apply the Extended 
Low Frequency Approximation, in our analysis we use a cutoff of 7 mHz as the trigger point for 
where we switched to the slower Rigid Adiabatic Approximation. This allows for quick and accurate 
frequency domain modeling for all but ~100 sources in our galactic simulations. 
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